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Abstract 

We present an algorithm for the adaptive tetrahedral integration over the Brillouin zone of crystalline materials, and apply it to 
compute the optical conductivity, dc conductivity, and thermopower. For these quantities, whose contributions are often localized in 
small portions of the Brillouin zone, adaptive integration is especially relevant. Our implementation, the woptic package, is tied 
into the wien2wannier framework and allows including a many-body self energy, e.g. from dynamical mean-field theory (DMFT). 
Wannier functions and dipole matrix elements are computed with the DFT package WiEN2k and Wannier90. For illustration, we show 
DFT results for fcc-Al and DMFT results for the correlated metal SrV 03 . 
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1. Introduction 

The theoretical description of crystalline solids is greatly 
simplified by their periodicity. The Bloch theorem for non-inter- 
acting electrons allows one to replace a sum over infinitely many 
discrete lattice vectors by an integral over a continuous k-vector 
restricted to the Brillouin zone (BZ). A similar simplification is 
possible in interacting systems where the crystal momentum k 
is a conserved quantity for a variety of excitations. This makes 
BZ-integration an indispensable part of any numerical technique 
for periodic solids. To evaluate the integral numerically, we must 
discretize the BZ in some way. The usual methods used in band- 
structure calculations rely on an a priori choice of the k-mesh 
which covers the BZ uniformly; e.g. using a straightforward 
summation [1-3], or the more sophisticated tetrahedron method 
[4]. This is a natural choice for the calculation of quantities 
such as the charge density, to which all k-points contribute. On 
the other hand, the transport or low-energy spectral properties 
are usually dominated by certain regions of the BZ, e.g. the 
vicinity of the Fermi surface. To compute such quantities, an 
inhomogeneous k-mesh adapted to a particular material may be 
a better choice. 

In the present article, we describe a technique to recur¬ 
sively generate an inhomogeneous k-mesh for periodic solids 
in three dimensions. Our implementation, the woptic package, 
is designed to calculate the optical conductivity, dc conduc¬ 
tivity, and thermopower of interacting electrons. However, the 
adaptive k-mesh management is encapsulated in a subprogram 
(ref ine_tetra) which may easily be adapted to other quanti¬ 
ties. 


Woptic operates in the context of dynamical mean-field the¬ 
ory (DMFT) [5, 6] for real materials. This “DFT-l-DMFT” ap¬ 
proach [7, 8] uses the band structure from density-functional 
theory (DFT) in the local-density or generalized gradient approx¬ 
imations (GGA) to construct an effective multiband Hubbard 
model, which is analyzed using the DMFT technique.^ Calcula¬ 
tion of the optical conductivity represents a post-processing step 
in this scheme. In the present implementation, we use inputs gen¬ 
erated by the WiEN2k [10], wien2wannier [11], and Wannier90 
[12] codes and a self energy (on the real-u; axis) from any DMFT 
solver. So far, only local self energies ^(u;) are implemented, but 
the approach allows including any self energy on top of WiEN2k. 
In particular, extension to a non-local S(A:, of) (e.g. from GW 
[13, 14], or from extensions of DMFT [15-17]) is simple as long 
as E(A:, (jj) may be obtained at any k. 

When vertex corrections are neglected (they are strongly sup¬ 
pressed in DMFT, see Sec. 2), the optical conductivity involves 
a BZ sum of contributions obtained from the k-resolved one- 
particle spectral functions and dipole matrix elements between 
the corresponding wave functions. We start by evaluating the 
optical conductivity on a uniform tetrahedral k-mesh. Next, we 
offer a refinement, test the sensitivity of the studied quantity and 
decide, for each tetrahedron, whether the refinement should be 
accepted. Accepting a refinement leads to the recursive genera¬ 
tion of additional k-points in a way that ensures the numerical 
stability of the algorithm. 

While the band structure part of the calculation uses aug¬ 
mented plane waves, the Hubbard model is naturally formu¬ 
lated in terms of localized orbitals. The transformation between 


- ^The calculations reported in this article use the GGA in the form of the PBE 

^Corresponding author functional [9]. 
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Figure 1: Schematic work flow of the woptic algorithm with an adaptive tetra¬ 
hedral mesh for the Brillouin zone integration. The algorithm is implemented 
in two main programs wopticjnain (evaluates the optical conductivity for a 
given k-mesh) and ref ine_tetra (adaptively refines the k-mesh) which are 
called by the driver script woptic, together with several support programs. Wan- 
nier90 provides the real space hopping matrix H(R), whereas the dipole matrix 
elements are computed either by the WiEN2k programs lapwl and optic, or 
by the woptic program convert_vr. See woptic user’s guide for more detailed 
work flow diagrams. 


the two bases is accomplished by the Wannier construction 
[11,18, 19]. 

An overview of the work flow of the adaptive refinement 
program called woptic can be found in Fig. 1. The work flow 
can be summarized as 

0. Choose an initial k-mesh and set iteration ^ = 0 (see 
Sec. 3.1 for the formal definition of the mesh). 

1. Compute the dipole matrix elements v^^\k) of Eq. (2) (see 
Sec. 4 for details). 

2. Compute the optical conductivity for a given k-mesh 
(see Sec. 2 for the formula used to obtain cr^^^). Extract the 
information which regions of k-space have a large contri¬ 
bution to the integration error (see Sec. 3.2 for details on 
the numerical quadrature and how the error is estimated). 

3. Stop if the change of the optical conductivity with respect 

to the previous iteration is below a given tolerance,^ 
otherwise refine the k-mesh where necessary, thus obtain¬ 
ing a new k-mesh, and return to step 1 with ^ ^ + 1 (see 

Sec. 3.1 for information on the refinement process). 

A more detailed version of the work flow, with a description of 
available modes, can be found in Sec. 4.2. 

This paper is structured as follows: Eirst, in Sec. 2, we give 
a description of the specific numerical problem to compute the 
optical conductivity. In Sec. 3, we specify the tetrahedral mesh 
and the refinement strategy. Eurthermore, we survey the esti¬ 
mation of the integration error necessary to mark tetrahedra for 
refinement in Sec. 3.2 and depict methods to increase the numer¬ 
ical performance in Sec. 3.3. In Sec. 4, we focus on practical 


^In the current version of the code, convergence has to be checked manually. 


considerations such as the available modes in the program and 
more details of the work flow and show numerical tests. Einally, 
in Sec. 5 we present two applications, elementary aluminum and 
the vanadate SrV 03 . 


2. Problem statement 


The Kohn-Sham Hamiltonian H of DET is diagonal in the 
Bloch-wave basis. Hence, the corresponding optical conductivity 
(T can be written in terms of the dipole matrix elements and 6- 
functions [20], setting 


(Tap{0) = - 


BZ 


S(Sc(k) - S,(k) - Q) 


Q 




,(k)4Ak)- ( 1 ) 


Here, Q is the external frequency, e is the electronic charge, the 
sum is over conduction (c) and valence (v) electrons with energy 
£c/vk, while 

Kmi^) = -- {>l'nk\da\<lfmk} ( 2 ) 

tUq 

are the dipole matrix elements with the Bloch states if/nu, and 
me denotes the electron mass. In general, the resulting value for 
(T depends on the number of k-points used in the integration as 
well as the applied quadrature rule. However, the evaluation of 
the integrand in Eq. (1) is numerically cheap and one usually 
proceeds to uniformly refine a given k-mesh until convergence. 

In the following we assume that we have constructed a map¬ 
ping U{k) from N Bloch states computed by WiEN2k [21] to N 
maximally localized Wannier functions (WEs) with Wannier90 
[12]; thus we have the Hamiltonian in Wannier space 


H{k) = U\k)E{k)U{k) with EUk) = SnmSnk- 0) 


This means that the electronic structure is not described in 
terms of bands Snu^ ^ ^ ^, but by a Hamiltonian ma¬ 
trix H{k) € This is appropriate, for example, when a 

set of WEs is required as an input for many-body calculations 
such as DET-bDMET. Also in light of possible DET-bDMET 
applications of the algorithm, we will not consider the effects of 
vertex corrections. In fact, while they can significantly affect the 
DMET results for other response functions like the spin/charge 
susceptibilities [22, 23], their contribution to the DMET optical 
conductivity is strongly suppressed (and vanishes exactly in the 
single-band case) due to the k-structure of the electronic current 
operator [24, 25]. On the basis of these considerations, the fol¬ 
lowing general expression for the optical conductivity [25, 26] 
can be written 


O-aA^) = - 


BZ 


doj - 


Tr [v^^{k)A{k, oj + Q)v^f^{k)A{k, co )]. (4) 

Here, / is the Eermi function at a given temperature. 


v7(A:) = U^k) vl(A:) U^k) = - {wMwsk) (5) 
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the dipole matrix rotated to the basis of k-space WFs and 

A the matrix spectral function 

Amnik, ^ ( 6 ) 

defined via the Green’s function 

G{k, oj) = [oj- H{k) - E(jt, w)]”‘ (7) 

and the corresponding electronic self energy 2 (A:, oj), all of which 
are also taken to be in the Wannier basis. Whereas in DFT, 
Eq. (1) scales linearly with the number of included bands, the 
general version, i.e., Eq. (4), scales quadratically with the num¬ 
ber of orbitals. Thus, for a basis set in which H(k) is not diagonal, 
such as WFs, it is useful to reduce the number of k-space evalu¬ 
ations, while keeping the level of accuracy as close to the DFT 
formalism as possible. 

In the specific case of DFT-hDMFT calculations, the (k, oj)- 
resolution required to resolve all features of the integrand of 
Eq. (4) depends in particular on the imaginary part of the self 
energy 2(A:, cl>) = contained in the DMFT Green’s function 
G(k, oj). For small 3E(6t;), a fine (k, 6 t;)-mesh is needed, since 

Tr [v'^CA:) A(k, coQ) v^(k) A(k, co)] ( 8 ) 

becomes sharply peaked. In many systems, the values of 3E(6g)) 
vary significantly between the different orbital manifolds, in par¬ 
ticular when considering transitions between localized orbitals 
(which have been treated, e.g., with DMFT) and itinerant orbitals 
(whose description usually remains at the DFT or at the Hartree 
approximation level). Hence, a useful approach is to adapt the 
k-mesh to the problem under consideration and use a finer reso¬ 
lution only where it provides a substantial increase of accuracy. 
Though the present implementation assumes local self energies 
2(6g)), generalization to a k-dependent oj) is straightforward 
as long as S(A:, oj) can be obtained for any k-point in reciprocal 
space. 

3. Algorithmic details 

3.1. Tetrahedral mesh 

The non-uniform triangulation of a three-dimensional do¬ 
main represents a formidable numerical task. The concepts sur¬ 
veyed in this section are not new but combine various well- 
known methods, and below, we will concentrate on the defini¬ 
tions necessary for the rest of this paper. For a more formal 
and complete introduction to tetrahedral triangulation see for 
example Ref. [27]. 

We denote the set of k-points of a certain tetrahedral triangu¬ 
lation by and the corresponding set of tetrahedra by T. Then, 
Nk = |7C| is the total number of k-points, and contains as 
elements or nodes the 3D coordinates 

nm = \kx ky , 1 < m < N^. (9) 

Furthermore, T stores a list of vertices 

Tm = [Vl V 2 V 3 V 4 ], 1 <m< Nt (10) 




Figure 2: An example of a 2D triangulation with a hanging node (a) shown in red 
and marked 1 (in this example we only discuss the nodes interior to the picture). 
Upon refinement of element A shown in (a), one arrives at the triangulation (b) 
with two additional hanging nodes 2, 3. This triangulation violates the regularity 
condition, since there are two hanging nodes 1, 2 on the same edge. To make 
(b) regular we must also refine element B; this is the mesh closure for this 
triangulation. 

where vi,..., V 4 e and Nt = \T\ is the total number of tetra¬ 
hedra (in practice, we store more information for each tetrahe¬ 
dron than just the vertices, see Sec. 3.2). Thus, the four vertices 
Vl,..., V 4 are nodes that define a tetrahedron e T and are 
themselves elements of 7C. We will use the term vertex only 
in connection to a specific tetrahedron, while a node is a gen¬ 
eral element of the set of k-points. A special type of node is a 
so-called 

hanging node: e is a hanging node if it lies on an edge 

of an element T without being a vertex of T. 

In the 2D visualization Fig. 2(a), node 1 is not a vertex of 
the central triangle B and hence a hanging node. Of course, a 
hanging node is a vertex of other triangles (here. A). 

We require 7“ to fulfill the following two conditions: 

regularity: No element T e T has an edge with more than 
one hanging node (see Fig. 2; otherwise, perform mesh 
closure, see below). 

shape stability: For all tetrahedra T e T, there is a prede¬ 
fined constant cs such that the radius of the circumscribed 
sphere rj satisfies r^/\T\ < cs where |r| is the volume of 
T. 

A large amount of the algorithmic effort in woptic is focused 
on keeping T shape stable and regular on refinement. If the ratio 
r\l\T\ becomes large this indicates a highly distorted tetrahedron 
(also called a degenerate element) and the numerical error of the 
integration rules may become large. A highly non-regular mesh, 
on the other hand, means that nearby regions of k-space are 
resolved very differently, see e.g. Fig. 2(b). This often leads to 
unstable convergence rates since some features of the integrand 
may not be fully resolved. 

In contrast to the 2D case of triangles, the refinement of a 
tetrahedron into 8 sub-tetrahedra of equal size is not unique and, 
in general, the resulting elements cannot all be similar to the 
original tetrahedron. In the following, we will depict Ong’s idea 
[28] for a refinement strategy where the shape of the element 
is at least confined to two classes, and shape stability is thus 
guaranteed. This strategy is based on the triangulation of the 
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parallelepiped defined by the three reciprocal unit vectors of 
the BZ into six tetrahedra of equal volume [Kuhn triangulation, 
Fig. 3(a)]. The resulting tetrahedra fall into two classes, in the 
following denoted by class 1 and 2. Specifically, it can be proven 
[28, 29] that for the refinements shown in Fig. 3, the resulting 
8 new elements of a tetrahedron will again belong to one of the 
two classes. Of the 8 new elements, 4 will share a vertex with the 
original tetrahedron and the other 4 will form a central octahe¬ 
dron. The difference between the two rehnement methods shown 
in Fig. 3 is how the central octahedron is split into tetrahedra. 
Since the triangulation of the central octahedron is not unique 
there are other strategies to ensure shape stability, e.g. based on 
the numbering of the vertices [27] . In practice, ensuring shape 
stability of the mesh T amounts to book-keeping of the classes 
of the tetrahedra upon refinement. 

In order to satisfy the regularity of the mesh T upon refine¬ 
ment, we add an additional step after the standard refinement of 
elements: mesh closure [30]. This procedure refines the elements 
with edges where regularity is violated (see Fig. 2(b) for a 2D 
example of a case where mesh closure is required). This leads 
to neighboring tetrahedra that will only differ by one level of 
refinement, since any tetrahedron which is refined twice while 
all of its neighbors remain in the initial state will automatically 
produce two hanging nodes on an edge. In this case the second 
refinement would trigger a refinement of the neighboring tetrahe¬ 
dra as well. Thus, when moving through k-space, the refinement 
level changes “smoothly”, i.e., regions with very fine resolution 
will not adjoin regions with very coarse resolution. These ad¬ 
ditional refinements lead to a higher number of total tetrahedra 
with respect to runs without mesh closure. Experience shows 
however that the price paid in performance is acceptable, since 
mesh closure helps avoid unstable runs of the algorithm. 

Summarizing, we obtain the following refinement strategy, 
assuming we have a regular mesh from the ^-th iteration of 
woptic (see Fig. 1) and a list of tetrahedra marked for refinement 
(see next section). 

1. Refine all marked elements of according to Fig. 3(c) 
and (e) for class 1 tetrahedra, and according to Fig. 3(b) 
and (d) for class 2 tetrahedra to obtain a rehned mesh 

2. Scan for hanging nodes. 

3. Mark all elements of for refinement whose edges 
violate the regularity condition, i.e., have multiple hanging 
nodes. 

4. If there are marked elements return to 1 with 
otherwise continue with 

This procedure gives a tree-like nested algorithm of refinement 
steps, which allows us to store the hanging nodes of and 
pass them on to the next iteration of woptic. In woptic, the rehne¬ 
ment strategy described above is implemented by the program 
ref ine_tetra (see Fig. 1). 

3.2. Integration error estimation and refinement 

Having outlined the salient points of mesh management, 
let us now discuss how elements of a mesh T are chosen for 
rehnement in the hrst place. Since we aim at minimizing the 


(a) 




Figure 3: (a) Kuhn’s refinement of a parallelepiped into eight tetrahedra as used 
by the woptic program. Class 1 tetrahedra are marked by beige vertices while 
class 2 tetrahedra are marked in blue. Note that every element of class 1 is 
mapped onto an element of class 2 and vice versa if mirrored along the main 
diagonal plane of the cube. The tetrahedra are first refined in 4 tetrahedra similar 
to the original one (having therefore the same class) and in a central octahedron 
denoted with black vertices, see panel (b) and (c). Depending on the class, the 
central octahedron is further split into 4 tetrahedra where 2 elements fall into 
the same class as the original tetrahedron and 2 in the respective other class, see 
panel (d) and (e). 


numerical integration error, we have to estimate the error et 
which an element T contributes to the overall error 


^tot — 


f dkgik)- YjSt 

BZ 


( 11 ) 


where g{k) denotes the integrand of the optical conductivity 


BZ 


/(tu + Q)-/(tu) 
doj - 

Q 


Tr \v^{k)A{k, oj + Q)v^{k)A{k, oj)] 

=: J" dkg{k) 

BZ 


( 12 ) 


and gT denotes an adequate tetrahedral quadrature rule (in g{k) 
we omitted the dependence on the external frequency for 
the moment, see below). For the integration over the internal 
frequency oj, we use a straightforward summation, exploiting 
only the weight factor f{oj -f f2) - f(oj) to limit the range of 
integration. For the k-integration we use two different rules: 
First, a linear 4-point rule 

S^T (13) 

i=l 
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with vi,..., V 4 being the vertices of T. Second, rule (13) can 
also be applied to the refined elements as 

sf = (14) 

i=\ i=\ 

where vji is the i-th vertex of the tetrahedron Tj (j = 1 ,..., 8 ) 
obtained from a refinement of T as introduced in the previous 
section. In the implementation, we evaluate the function g(k) 
on the 10 points required to apply both rules Eq. (13) and (14). 
The 4 vertices of T required for the rule 4p are also included in 
Ti,..., Tg. The latter additionally have 6 midpoints on the edges 
of T. Thus, G 7“ is represented by the 10 nodes plus the class 
of the tetrahedron (which is required for the refinement strategy, 
see previous section). 


Tm = [vi V 2 V 3 V 4 ; nu ni 3 ^23 ^14 n23 ^ 34 ; 1 or 2], (15) 

where nij is the midpoint between the vertices v/ and v^ . Note 
that the nested nature of our quadrature rules allows us to re-use 
values of g(k) in following iterations of the algorithm. 

To estimate the contribution which T adds to the total er¬ 
ror Etot, we compare the results of Eqs. (13) and (14), which 
means that the same rule is compared for two different levels of 
refinement [31]. Thus, our error estimator is 


Et - 



4p 

\St 


4pr 

■ 


(16) 


Since the rule (14) is obtained by applying the rule (13) to the 
sub-elements Ti,..., Tg which would be new elements if T was 
refined, the error estimate et provides a measure of how much a 
refinement of T would improve the numerical integration. 

The dependence of the optical conductivity (Tapm) on the 
external frequency Q and the directional dependence {ap) have 
been neglected so far. To take these dependencies into account, 
all error estimates of an element are averaged, 

er = ^ T (O), ap € [xx, xy, xz, yy, yz, zz]. (17) 


To mark certain elements for refinement, we apply a standard 
procedure for adaptive mesh algorithms [30]: an element T is 
marked if 


> 0 maxGj’/, (18) 

T'eT 

where 0 G [0,1] is a parameter determining the harshness of 
the refinement. A value of 0 = 0 means that all elements sat¬ 
isfy (18), i.e. uniform refinement, whereas large values of 0 
lead to highly adaptive meshes. In woptic, the error estimation is 
partly performed by woptic jnain and partly by ref ine_tetra 
(see Fig. 1). The former computes the integrand, the latter calcu¬ 
lates the error estimators et and marks the elements for refine¬ 
ment according to Eq. (18). 

In metallic cases, the optical conductivity cr(Q) for ^ 0 
has a Drude contribution corresponding to a Lorentzian at Q = 0 


which is broadened by 3S(0). Thus, if 3E(0) is small, the error 
estimator (17) is often dominated by the values around the Fermi 
level Q = 0 and the algorithm mainly resolves the Fermi surface. 
This behavior may be adequate when one is interested in the dc- 
conductivity or the thermopower, but for the optical inter-orbital 
transitions at higher energies one might favor a better description 
in that region. In this case, another error estimator instead of 
Eq. (17) is more appropriate: 


1 


Q-ajS 


ap G {xx, xy, xz, yy, yz, zz}, (19) 


where the additional factor Q, attributes a larger weight to the 
error at higher frequencies. 


3.3. Performance and symmetry considerations 

Given that one evaluation of the function g{k) from Eq. (12) 
is numerically expensive, the number of total evaluations should 
be kept as small as possible. For this reason, we use two tech¬ 
niques: (/) re-using the data from previous iterations and (//) 
taking into account the symmetries of the crystal. The first point 
is the main reason for choosing the two nested quadrature rules 
Eqs. (13) and (14), since, as mentioned above, a refinement of an 
element T eT yields at most 6 new nodes. Moreover, neighbor¬ 
ing elements with similar refinement level share nodes with T. 
Thus, though our integration rules are of low order, they repre¬ 
sent an efficient choice in terms of the number of total evaluation 
points. 

To understand (//), i.e. how to increase the performance by 
symmetry, let us define matrices S c describing the sym¬ 
metry operations of the crystal in a Cartesian coordinate system. 
Furthermore, c denotes the symmetrized k-mesh, i.e. the 
reduced mesh when the symmetry operations of S are exploited, 
with the corresponding mapping ms '. n ns such that ^ g 
and ns g If one replaces each vertex v of each element 
of T by its reduced vertex m^(v), one formally obtains a new 
tetrahedral mesh Ts. Note that Ts might include elements that 
do not correspond to real tetrahedra but have e.g. equal nodes 
when multiple vertices of an element of T have been mapped 
onto the same k-point in the reduced set After the mapping 
T ^ Ts there are in general multiple occurrences of an element 
Ts. For simplicity, Ts in the following denotes the reduced sym¬ 
metrized mesh, i.e. all elements Ts e Ts are only considered 
once and carry a weight w , which accounts for the volume and 
multiplicity of Ts. 

The numerical quadrature to yield the optical conductivity 
according to Eq. (11) is given by 

<T=y^gT. ( 20 ) 

TeT 

It is important to stress here that one cannot simply replace 
n e T within the rule gr by ms(n) G since cr and gr are 
tensors. Hence, rotated quantities have to be used: 

^ Z Z (21) 

' ^ W,eST,eTs 


5 





Hw{k) 

E’>{k) 


formation from the initial, “random” gauge^ to the new gauge 
[19]. One then expresses the Kohn-Sham Hamiltonian and the 
dipole matrix elements in the new gauge using Eqs. (3) and (5), 
to repeat: 

H(k) = U^(k) E(k) U(k) with Enm(k) = Snm^nk, 
v^(k) = U^(k)v(k) U(k). 

The smoothness of the the Wannier gauge guarantees that the 
corresponding basis functions (i.e. the WFs) are exponentially 
localized, and hence that the Fourier transforms of the Hamilto¬ 
nian 


Figure 4: Possible optical transitions and corresponding Hamiltonians according 
to Eq. (26) for various manifolds in a low energy model for SrV 03 . In the 
schematic visualization of the spectrum, the t 2 g manifold where the Wannier 
projection is performed is marked Hw(k) and other parts of the spectrum 


In this approach it is sufficient to compute g(ns) for Us e and 
this, depending on the symmetry of the problem, may yield a 
considerable speed-up. 


HrsiR) = 


E H(k)e = {w,i,\H\WsR) 


( 22 ) 


and dipole matrix elements 


V*(^) = E 


- - {Wro\da\WsR} 

trie 


(23) 


4. Practical usage 

After the general description of the algorithm in the previous 
section, let us now turn to the connections of woptic to the 
program packages WiEN2k and Wannier90, the available modes 
of operation and a more detailed work flow. As a prerequisite 
to start a woptic calculation, two other packages are needed: 
WiEN2k [21] (including wien2wannier [11]) and Wannier90 [12]. 
This set of programs allows constructing maximally localized 
WFs from WiEN2k (see Refs. [11, 32] for a detailed description). 

To employ the adaptive integration described in the previ¬ 
ous section, we have to be able to generate the dipole matrix 
v'^(k) and the Hamiltonian H(k) at arbitrary, a priori unknown, 
k-points. Two approaches are implemented in woptic: (1) matrix 
elements between WFs can be obtained by Fourier transforma¬ 
tion from direct space — this is called called interp mode in 
the program; or (2) these matrix elements can be recomputed by 
WiEN2k every iteration — optic mode. 

Two issues are central in the following: 

Gauge invariance of the optical conductivity (4) and other ob¬ 
servables, i.e. the independence of the random gauge of 
the Bloch states. While the trace in Fq. (4) is gauge in¬ 
variant, its building blocks v'^(k) and A{k, oS) are not. It 
is therefore essential that they be expressed in the same 
gauge. 

Mixed transitions, dipole transitions between Wannier states 
on the one hand, and Bloch states that are not included in 
the initial Wannier projection on the other. 

Gauge invariance is easily ensured if approach 1 (interp 
mode) is available. The Wannier construction consists in finding 
a k-smooth gauge on the initial k-mesh TCw This amounts to 
finding a set of unitary matrices U{k) which represent the trans- 


converge rapidly with the mesh spacing of TCw (Here and below, 
Wruir) is a direct-space WF.) 

Thus we obtain direct space hopping and dipole parameters 
for all the significant neighbor shells (exponentially small contri¬ 
butions of the tails of the WFs are neglected). Having the direct 
space representation allows us to compute the matrix elements 
at an arbitrary k-point q outside TCw via 

H(9)= y] H(/?)e‘«A (24) 

V*(^) = ^ v*(/?)e‘«*. (25) 

To emphasize, the R-sums run over the set of lattice vectors 
dual to TCw, but the localization of the WFs allows us to apply the 
reverse Fourier transform for q ^ TCw with excellent accuracy. 
This procedure is known as Wannier interpolation [33, 34]. 

Mixed transitions arise when we want to extend the optical 
conductivity, which we have written in Fq. (4) in terms of a 
trace over the WFs (inner or Wannier window), to include states 
not covered by the Wannier projection (outer or Bloch window). 
The motivation for the outer window is normally to compute the 
optical conductivity over a larger frequency range (see Fig. 4 for 
a visualization of the model and the possible transitions). 

For a DFT-hDMFT calculation, the self energy E(6 g)), describ¬ 
ing correlation effects beyond GGA, has to be provided on the 
real axis. In the following, this situation is referred to as inter¬ 
acting. Woptic also provides the option to set the self energy to 
a small imaginary constant 1^(oj) = -iS to mimic broadening, 
e.g. from impurity scattering."^ We will refer to bands treated 


^More precisely, the gauge determined by diagonalization of H(k) in the 
electronic structure code. 

^This means that in the non-interacting case, the lifetime broadening (result¬ 
ing also in a finite width of the Drude peak) is added “by hand”, while in the 
interacting case, it arises naturally from DMFT. 
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in this way as non-interacting. By supposition, the Bloch states 
in the outer window are adequately described in DFT and thus 
non-interacting in this sense. 

When an outer window is included, we split the Hamiltonian, 
Wannier transformation, and dipole matrices into Wannier, Bloch, 
and mixed parts. Likewise, the trace in the optical conductivity 
(4) splits into Wannier-Wannier, Wannier-Bloch, and Bloch- 
Bloch terms. Denoting the diagonal energy matrix connected to 
the states above (below) the WFs in energy by we can 

write the “large” Hamiltonian, which is block-diagonal. 


(E\k) 


mk) = 


H(k) 


E%k) 


(26) 


Analogously, the large Wannier transformation matrix is 


(1 


mk) = 


U(k) 

V 


(27) 


It affects only the inner window and leaves the outer window un¬ 
changed. Additionally, we have the large dipole matrix "V^^ik) = 
{^mk\daVknk) with the indices n,m running now over all 
bands in the outer window instead of only the inner window. 
Inserting E((k) into Eqs. (7) and (6) yields the block-diagonal 
matrix spectral function ^(k, oj) of the large system. Together 
with the large dipole matrix in the Wannier basis 


^'^{k) = V^{k)^{k)V{k) 



'(vp 




(v“) 

(V*) 




(Ki) 



(28) 


this spectral function can be used in Eq. (4) to give a more 
complete description of optical transitions in the system. 

Approach 1 (interp) is straightforward for the Wannier- 
Wannier transitions, but not directly applicable to the mixed 
dipole matrix elements 


= L {Wrk\daWik) (29) 

n 


because their Fourier transform v^.(/?) does not decay with \R\. 
To salvage Wannier interpolation in the presence of an outer 
window, we define the quantity 

w“f (*, = L (30) 


where the index i runs over all non-Wannier states (which are 
non-interacting, hence their matrix spectral function is di¬ 
agonal). Its Fourier transform w^^(R, oj) decays with \R\, albeit 
more slowly than v'^(R) [35]. 

With these two interpolated quantities, and using the WiEN2k 
programs lapwl and optic [20] to compute the Bloch energies 
E^^^(q) and Bloch-Bloch dipole matrix elements v^j(q), we can 


evaluate the trace (8) at any new k-point. On the other hand, 
the interpolation errors from w^^(q, oj) may get large, see next 
section and Ref. [35] for tests. (Interpolation errors from v'^(^) 
are insignificant so long as properly localized WFs are found.) 

As an alternative in the mixed case, we turn to approach 2 
(optic mode): computing "Viq) and E^^^(q) at new k-points us¬ 
ing lapwl and optic. (The Hamiltonian H(q) is still computed 
using Eq. (24).) Because the (inner-window) spectral function 
A(q, oj) is computed in the Wannier basis, but optic yields the 
dipole matrix elements in the Bloch basis, we need the Wannier 
transformation U(q) on the new k-points to mediate between 
them. Since the Bloch basis is the one in which the Hamiltonian 
is diagonal, we can obtain U(q) by diagonalizing H(k) (inverting 
Eq. (3)). 

The problem with approach 2 (optic) lies in the arbitrari¬ 
ness of the Bloch gauge. Since the U(q) obtained via Wannier 
interpolation are computed by diagonalization, they are deter¬ 
mined only up to the phases of the respective eigenvectors. In 
the non-interacting case, where the matrix spectral functions are 
diagonal, these phases evidently cancel in the trace (8); in fact, 
this reasoning can be extended to the interacting case as long as 
the Wannier self energy is a scalar (diagonal in and independent 
of the orbital index), e.g., because of crystal symmetry. 

From a different point of view, the Bloch states if/nk, obtained 
as solutions of independent eigenproblems at each k-point, carry 
“random” phases. The original U{k) from Wannier90 take these 
phases into account in constructing smooth functions Wrkir) = 
Y,n k^nr(k) ^nki^) of k. Thcsc phascs are included both in U{k) 
and v^{k) and hence cancel when calculating the dipole matrix 
elements in Wannier space using Eq. (28). 

However, if the adaptive k-mesh algorithm now selects a new 
point q, the phase of if/nq is included in the recalculated v^{q) but 
not in U{q) obtained as explained above. Hence, the “random” 
phase may enter into the trace (8) both through the Wannier- 
Wannier and the mixed transitions (for the Bloch-Bloch transi¬ 
tions it cancels). The resulting random-gauge problem leads to 
errors in the results whose magnitude is a priori unknown. 

So far we have assumed that the transformation between 
the Bloch and Wannier states at each k-point is accomplished 
by a unitary matrix U{k). This excludes the disentanglement 
procedure [36] implemented in Wannier90, where additionally 
a rectangular matrix V(k) intervenes. In fact, woptic supports 
disentanglement only in interp mode without an outer win¬ 
dow (Wannier-Wannier transitions only). It is not clear how the 
method may be extended to the general disentangled case [35]. 

4.1. Benchmarks of interpolation and random-gauge errors 

In Fig. 5, we compare the optical conductivity of SrV 03 
from the optic and interp modes, including a self energy 
from DMFT on the V-t 2 ^ states (see Sec. 5.2 for details on 
the DMFT calculation). Since in this material the self energy 
is orbital independent by symmetry, there is no random-gauge 
problem, and this case can be regarded as a test for the Wannier 
interpolation of and 

In order to quantify the random-gauge errors, we require a 
test case which is complementary to the one of Fig. 5, i.e., where 
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Figure 5: Optical conductivity of SrV 03 computed in the optic (red curve) 
and interp (blue curve) modes. For this model, whose self energy is orbital- 
independent due to the cubic symmetry, optic mode can be considered exact. 
Thus the difference between the curves reflects the interpolation errors afflicting 
the mixed (Wannier-Bloch) transitions in interp mode. Note that the outer 
window of included bands is smaller here than in Fig. 9, which is why the optical 
conductivity starts to drop off above cd ^ 7 eV. Here, the Wannier projection 
comprises the 3 V-tig bands and 14 bands are included in the outer window, 
corresponding to O-p and V-d states. 

interpolation is reliable, but the gauge problem is in effect. To 
this end, we construct a Wannier projection covering exactly 
those bands of SrV 03 included in Fig. 5, and apply the ^-hg 
self-energy from Fig. 5 to the t 2 ^-derived orbitals of this larger 
projection.^ Since all included transitions are now between WFs, 
we need to interpolate only and can rely on interp mode 
as a reference; but since we apply a non-trivial self energy only 
to the t 2 ^ orbitals, is now strongly orbital-dependent and 
optic mode suffers from the gauge problem. The results are 
shown in Fig. 6. 

Comparing Figs. 5 and 6, we find that the errors from 
interpolation and from the random-gauge problem are similar in 
magnitude. In both cases, the qualitative features are preserved. 
The quantitative differences must be viewed in relation to other 
sources of uncertainty in these calculations. Most importantly, 
in DMFT the self energy on the real-u; axis is typically obtained 
through analytic continuation from the imaginary-u; axis. This 
leads to uncertainties which can easily be comparable to the 
errors observed in Figs. 5 and 6. It also bears mentioning that the 
orbital symmetry which protects optic mode from the random- 
gauge problem is broken especially sharply in the 14-band model 
used above (nontrivial on the t 2 ^ orbitals, S/ = -id on the 
others). 

Fig. 6 contains a third curve, which corresponds to interp 
mode on a model constructed with disentanglement. At the R 
point of the BZ, the V-e^ bands are entangled with Sr-s bands, 
as seen in the band structure of Fig. 4. These unwanted bands 
can be removed using disentanglement [36]. The corresponding 
curve in Fig. 6 is practically identical to the one without disen¬ 
tanglement, except for the region around 7 eV, where transitions 
involving the entangled states are relevant. 


^This is unphysical but yields a convenient test case. 



Figure 6: Optical conductivity of a model derived from SrV 03 by imposing the 
self energy from the 3-band model on the t 2 g-like orbitals in a 14-band Wannier 
projection, which includes the same bands as in Fig. 5: O-p, V-t 2 g, and V-e^. For 
this model without mixed transitions, interp mode can be considered exact, thus 
the difference reflects the random-gauge errors in optic mode. The dashed blue 
curve is from interp mode using a Wannier projection with disentanglement of 
the Sr-s bands crossing the V-e^ bands at the R point of the Brillouin zone [35]. 
Since this corresponds to a somewhat different model, complete agreement with 
the solid blue curve cannot be expected. 

For details on the issues discussed in the preceding para¬ 
graphs (to wit: the random-gauge problem in optic mode, in¬ 
terpolation of and and disentanglement in relation to 
woptic), we refer to Ref. [35]. 

4.2. Program details 

To summarize, woptic offers two main choices for the ma¬ 
trix element mode (corresponding to the program parameter 
matelmode). This choice determines the method to compute 
the dipole matrix elements v{k). The modes and corresponding 
keywords in the woptic input file are: 

Wannier interpolated dipole matrix elements (interp): 

Apply Wannier interpolation (25) to the dipole matrix 
elements v^^{k) in the Wannier gauge (5), as well as to 
the Hamiltonian (3). In the presence of an outer window, 
the mixed transitions are calculated from the quantity 
Wrsik, oj) (30), which is likewise interpolated [35]. 

Ab initio dipole matrix elements (optic): Obtain the dipole 
matrix elements vf.(^) in the Bloch gauge (2) from the 
WiEN2k programs lapwl and optic; the Hamiltonian 
H(q) from Wannier interpolation; and the transformation 
U(q) to the Wannier gauge by diagonalization of H(q). 

Interp mode is reliable for the Wannier-Wannier and Bloch- 
Bloch transitions, but the Wannier-Bloch terms acquire interpo¬ 
lation errors. Optic mode is reliable whenever the self-energy in 
the Wannier gauge is diagonal and orbital independent (by sym¬ 
metry, or in the non-interacting case), but the Wannier-Wannier 
and Wannier-Bloch terms acquire errors due to the random- 
gauge problem when it is not. In our tests, the errors from these 
two issues are comparable. 

We conclude this section with a summary of the detailed 
work flow of woptic. 
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0. A preliminary run of ref ine_tetra prepares the initial k- 
mesh and the initial set of tetrahedra In interp 
mode, obtain v^{R) and w^^(R, oj) (if mixed transitions 
are requested) from Eqs. (5) & (30). Set ^ = 0. 

1. Determine which k-points of were not in i.e., 

find \ 

2. Obtain the Hamiltonian H(k) for all k-points in \ 

yia Eq. (24). 

3. Obtain the dipole matrix v^nik) for all k-points in \ 

from the WiEN2k programs lapwl and optic in 
optic mode, or from Eq. (25) in interp mode. 

4. Call wopticjnain 

(a) Rotate v{k) to the Wannier basis (Eq. (5)). 

(b) Load the self energy 2(6t;) and determine the Green’s 
function G{k, oj) according to Eq. (7) for all k-points 
in \ 

(c) Evaluate the contributions to the optical conductivity 
g{k) from Eq. (12) for all k-points in \ 

and load the old data g{k) for 

(d) Perform tetrahedral integration for using Eq. 
(21) and obtain the optical conductivity cr^^^ipj) 

5. Call ref ine_tetra 

(a) Determine the error estimators for T e via 
Eqs. (13 & 14) and (16 & 17), respectively. 

(b) Mark the elements of for refinement if they 
satisfy the criterion (18), obtaining c 

(c) If is empty, i.e. no elements have been marked, 

set = 7“^^^ (and = 7C^^^) and exit from 

ref ine_tetra. 

(d) Refine the marked tetrahedra 7“^ according to their 
class and the rules shown in Eig. 3, leading to the 
refined mesh of and a new set of k-points 

^ref 

(e) Perform mesh closure: Mark the non-refined tetrahe¬ 
dra of for refinement if they violate the regular¬ 
ity condition, i.e. if they have more than_^one hanging 
node on an edge, and obtain . If is empty, 
the mesh is regular: set 7“^^+b = 

and exit. Otherwise set = Tln \ 
rj^{£) _ <^^(0 and return to step 5d. 

6. If ^ < ^niax return to step 1; otherwise, exit. 

The key difference between the modes is in step 3, where the 
matrix elements for the new k-points are computed. 

5. Applications 

5.7. Aluminum 

As a simple example to show that our adaptive procedure can 
reproduce standard WiEN2k results at much lower computational 
cost we have chosen fcc-Al. As shown before [20, 38], there is a 
strong dependency of the optical properties on the quality of the 
BZ integration. Using a regular grid and the tetrahedron method 
as implemented in WiEN2k’s optic module, up to 20000 k- 
points in the irreducible wedge of the BZ are necessary to obtain 
converged results for the optical conductivity. With our adaptive 


mesh the number of k-points can be greatly reduced. Note that 
(t{oj) also depends crucially on the applied broadening scheme, 
and small differences may appear between the two methods, 
because they are somewhat different in the way the broadening 
is introduced. The optic module in WiEN2k first calculates the 
imaginary part of the unbroadened dielectric function, 62 ( 0 ;), 
due to interband transitions and the plasma frequency due to 
intraband transitions, and adds smearing later (using Lorentzian 
broadening). On the other hand, the Green’s function method 
uses a related broadening constant 6 for the self energy ^(u;) = 
-\d which enters into the Green’s function (7). 

The reason for the slow convergence of cr with the number 
of k-points is quite obvious when one inspects the band structure. 
Consider the first interband peak at 61 ; 1.5 eV (Eig. 7(a)). In 

the relevant energy range, only a narrow region of the BZ around 
W contributes, as seen in the k-resolved contributions to (t{oj) 
(Eig. 7(b)). Comparison with the bandstructure (Eig. 7(c)) shows 
that these contributions stem from four bands near the Eermi 
level at W (two particle and two hole bands). A regular mesh 
must be quite dense to properly sample these small portions of 
the BZ, while our adaptive scheme is much more efficient. This 
is illustrated in Fig. 8 , where one can see that the regions around 
the W-point are refined most and have the smallest tetrahedra. 

5.2. Strontium vanadate 

As a second application of the method, we present in this 
section calculations for SrV 03 . Its low-energy electronic struc¬ 
ture is dominated by the degenerate 3d-t2^ orbitals of vanadium 
and it constitutes a textbook example of a strongly correlated 
metal, perfectly suited for illustrating the intermediate steps and 
the final results of the woptic package. Specific details about the 
crystal and electronic structures of SrV 03 can be found, e.g., in 
Refs. [26, 32]. In fact, SrV 03 is a good example of a situation 
where DFT cannot accurately describe the low-energy optical 
response of the system,^ necessitating the inclusion of local cor¬ 
relation effects, e.g. by means of DMFT. In this work, we focus 
on the workings of the adaptive algorithm rather than physical 
implications of our results; see Ref. [26] for a discussion about 
the latter. 

In Fig. 9, the intermediate (^ = 1) and final {i = 5) results for 
the optical conductivity of SrV 03 are reported, in each case for 
an interacting and a non-interacting calculation. Experimental 
results from Ref. [39] are reproduced for comparison. Here, 
the DMFT calculations have been performed using a Kanamori 
interaction [42] with parameters U = 5.05 eV, 77' = 3.55 eV 
and J = 0.75 eV, consistent with the setup used in Ref. [26]. 
Since the random-gauge problem is absent in this case due to 
the cubic symmetry (see Sec. 4), we use optic mode for the 
interacting optical conductivity. One immediately notes the role 
of the electronic correlations, as evidenced by a significant shift 
of optical spectral weight from the Drude peak and the frequency 
window between 3 and 5 eV to higher energies as compared to 


^Another prototypical situation is the analysis of the optical spectroscopy 
experiments in V 2 O 3 [26, 32, 40, 41], where the corrections generated by the 
inclusion of electronic correlation in DFT+DMFT are even larger than for the 
SrV 03 case. 
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(a) 


(b) 



Figure 7: The non-interacting optical conductivity of A1 computed by woptic compared to corresponding results from the WiEN2k optic package and experimental 
data [37] in panel (a). For convergence of the uniform WiEN2k calculation a large number of (symmetrized) k-points is required, while due to the adaptivity, woptic 
converges for a much smaller number of k-points. After convergence both programs yield similar results, in particular the experimental peak position is reproduced. 
(An arbitrary scaling has been applied to the experimental results.) The contributions to the optical conductivity cr resolved in {k, (u)-space (b). Only a small part of 
k-space around W contributes to cr and the contributions can be understood by identifying possible transitions in the band structure (c). 
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c/} 
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^v-axis 


kjc-axis 


ky-axis 


kx-axis 




Figure 8: The unsymmetrized initial tetrahedral mesh used in the calculations for elementary A1 with 3072 tetrahedra and 4913 k-points (left). The highly 
adaptive mesh after 6 iterations 7“^^^ with 13152 tetrahedra and 24667 k-points (right). For the sake of visualization only the smallest tetrahedra are shown. The 
region yielding the largest contribution to the estimator (17) is located close to the W-point around [0.7500 0.5000 0.2812] in terms of the primitive vectors of the 
reciprocal unit cell. 
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Figure 9: The interacting and non-interacting optical conductivity of SrV 03 computed by the adaptive algorithm woptic (optic mode) compared to experiment [39] 
in panel (a). While the interacting result is much closer to experiment, the convergence of the algorithm requires far fewer iterations than for A1 both in the interacting 
and in the non-interacting case. This indicates that larger regions of k-space contribute to the optical conductivity, as also suggested by the contributions to cr resolved 
in (k, 6cj)-space (b). (Note also the different scales of the contributions here and in Fig. 10.) A significant part of the optical conductivity in the energy window under 
investigation stems from O-p ^V-d transitions, which can be understood by inspection of the band structure (c). In panels (a) and (b), the dashed line indicates the 
restricted frequency window (0 to 5 eV), which has been considered for the illustration of the evolution of the tetrahedral mesh in Fig. 10. In the non-interacting 
case, the broadening, in particular the finite width of the Drude peak, arises solely from the small imaginary self energy E = -id which is added for this purpose. By 
contrast, in the interacting case, it is a direct result of DMFT. 



^X-axis 


1 



^3,-axis 


^;c-axis 


Figure 10: The unsymmetrized initial tetrahedral mesh 7“^^^ used in the non-interacting calculations for SrV 03 with 3072 tetrahedra and 4913 k-points (left). After 5 
iterations, in the highly adaptive mesh regime (6 = 0.9), the mesh 7“^^^ has 4920 elements and 9533 k-points (right); for the sake of visualization only the smallest 
tetrahedra are shown) 
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the non-interacting calculations. Such many-body effects are 
expected due to the correlated nature of the 3d-t2^ orbitals of 
SrV 03 , and they significantly improve the match between theory 
and experiment. 

A more detailed analysis of the woptic results in Fig. 9 shows 
that the convergence of the adaptive algorithm is much faster than 
for the A1 case of the previous section. This can be understood by 
plotting the k-resolved contributions to cr(oj) (i.e., the integrand 
of the k-summation of Eq. 4), as shown in Fig. 9(b): In the wide 
energy range considered (i.e., up to 15 eV), the contributions 
to the optical conductivity are spread over all of k-space. As a 
consequence, the adaptivity of the woptic algorithm becomes 
less important and the final adaptive mesh (not shown) essentially 
coincides with one obtained in a uniform calculation. This would 
be different when focusing on the low-energy region (e.g., up to 
5 eV as marked by the dashed lines in Fig. 9) - as is usual when 
comparing with optical spectroscopic experiments. In that case, 
the predominant contribution to cr(oj) below 5 eV (apart from 
the Drude peak) is from around the T point, and between T and 
X. This corresponds to the peak in cr(oj) located at about 4 eV in 
the non-interacting spectrum, which originates not from optical 
transitions within the V-t 2 ^ orbitals, but rather from transitions 
between the 0-2p and the V-e^ bands. 

This situation is well reflected in the evolution of the tetrahe¬ 
dral mesh, reported in Fig. 10 for a calculation in the window 
up to 5 eV, and performed with a highly adaptive mesh (6 = 0.9) 
for the sake of illustration. In the right panel of Fig. 10 the re¬ 
sulting tetrahedral mesh after 5 iterations is visualized, showing 
refinements essentially from T ^ X, which is exactly the region 
we identified in Fig. 9(b) and (c). Moreover, especially along 
this k-path, the p and the e^ orbitals mainly responsible for the 
optical transitions are relatively flat and hence difficult to resolve 
in (k, 6c))-space, which explains woptic’s behavior in this case. 

6. Summary 

We have developed a flexible and efficient adaptive BZ inte¬ 
gration algorithm based on a recursively generated tetrahedral 
k-mesh. In regions where the numerical error would otherwise 
be large, the k-mesh becomes fine, whereas it remains coarser 
elsewhere. We apply this approach to the optical conductivity in 
a Wannier basis, with the possibility to include a many-body self 
energy on top of the DFT band structure. The peakedness 
of the contributions in k-space is determined mainly by the band 
structure and by the imaginary part of the self energy, which 
broadens the peaks. Thus, weakly interacting materials, where 

is small, tend to have more sharply peaked contributions. 

Results for A1 and SrV 03 illustrate the algorithm and its 
performance. These calculations would require much more com¬ 
putational effort using uniform k-grids. 

The woptic package, our implementation of the adaptive k- 
mesh algorithm in the framework of WiEN2k, Wannier90, and 
DMFT, is available at http: //woptic. github. io. In addi¬ 
tion to the ready-made computation of the optical conductivity, 
dc conductivity, and thermopower, the k-mesh management code 
may easily be adapted to other quantities, in particular where a 


conventional tetrahedron integration is impossible or impracti¬ 
cal. 

In order to include transitions involving bands beyond the 
Wannier projection, an outer band window may be defined, al¬ 
though this leads to certain numerical problems in some cases 
(see Sec. 4). The outer bands will be described at the DFT level, 
i.e. without a self energy. Apart from the physical, k-integrated 
quantities, woptic also provides tools to examine the k-dependent 
contributions (as in Figs. 7 and 9 (b)), which often provide valu¬ 
able physical insight. 
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